model { for(i in 1:N){ mu_M[i] <- alpha*x[i] M[i] ~ dnorm(mu_M[i],prec1) mu_y[i]<- c*x[i]+beta*M[i] y[i] ~ dnorm(mu_y[i],prec2) mu_y3[i]<-beta*M[i] #when x=0 mu_y4[i]<-c+beta*M[i] #when x=1 de[i]<-(mu_y4[i]-mu_y3[i])/deltax mu_y2[i]<-c*x[i]+beta*(M[i]+deltam) ie[i]<-alpha/deltax*(mu_y2[i]-mu_y[i])/deltam #alpha is alpha(1-0) te[i]<-ie[i]+de[i] } alpha ~ dnorm(0.0,0.000001) beta ~ dnorm(0.0,0.000001) c ~ dnorm(0.0,0.000001) var1 ~ dgamma(1,0.1) var2 ~ dgamma(1,0.1) prec1 <-1/var1 prec2 <-1/var2 }